Functions and Data Load
loadPkg = function(pkg, character.only = FALSE) {
if (!character.only) { pkg <- as.character(substitute(pkg)) }
if (!require(pkg,character.only=T, quietly =T)) { install.packages(pkg,dep=T,repos="http://cran.us.r-project.org"); if(!require(pkg,character.only=T)) stop("Package not found") }
}
loadPkg(knitr)# Fix outliers
outlierKD2 <- function(df, var, rm=FALSE) {
#' Original outlierKD functino by By Klodian Dhana,
#' https://www.r-bloggers.com/identify-describe-plot-and-remove-the-outliers-from-the-dataset/
#' Modified to have third argument for removing outliers inwtead of interactive prompt,
#' and after removing outlier, original df will not be changed. The function returns the new df,
#' which can be saved as original df name if desired.
#' Check outliers, and option to remove them, save as a new dataframe.
#' @param df The dataframe.
#' @param var The variable in the dataframe to be checked for outliers
#' @param rm Boolean. Whether to remove outliers or not.
#' @return The dataframe with outliers replaced by NA if rm==TRUE, or df if nothing changed
#' @examples
#' outlierKD2(mydf, height, FALSE)
#' mydf = outlierKD2(mydf, height, TRUE)
#' mydfnew = outlierKD2(mydf, height, TRUE)
dt = df # duplicate the dataframe for potential alteration
var_name <- eval(substitute(var),eval(dt))
na1 <- sum(is.na(var_name))
m1 <- mean(var_name, na.rm = T)
par(mfrow=c(2, 2), oma=c(0,0,3,0))
boxplot(var_name, main="With outliers")
hist(var_name, main="With outliers", xlab=NA, ylab=NA)
outlier <- boxplot.stats(var_name)$out
mo <- mean(outlier)
var_name <- ifelse(var_name %in% outlier, NA, var_name)
boxplot(var_name, main="Without outliers")
hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
title("Outlier Check", outer=TRUE)
na2 <- sum(is.na(var_name))
cat("Outliers identified:", na2 - na1, "\n")
cat("Propotion (%) of outliers:", round((na2 - na1) / sum(!is.na(var_name))*100, 1), "\n")
cat("Mean of the outliers:", round(mo, 2), "\n")
m2 <- mean(var_name, na.rm = T)
cat("Mean without removing outliers:", round(m1, 2), "\n")
cat("Mean if we remove outliers:", round(m2, 2), "\n")
# response <- readline(prompt="Do you want to remove outliers and to replace with NA? [yes/no]: ")
# if(response == "y" | response == "yes"){
if(rm){
dt[as.character(substitute(var))] <- invisible(var_name)
#assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
cat("Outliers successfully removed", "\n")
return(invisible(dt))
} else {
cat("Nothing changed", "\n")
return(invisible(df))
}
}loadPkg('dplyr')##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionloadPkg('readxl')
loadPkg('readr')
#Read in general census data
census_2015 <- read.csv(file = 'acs2015updated.csv')
#Read in Specifically Education Data
edu_data_2015 <- read_csv("ACS_15_5YR_B15003_renamed.csv", skip = 1) %>% mutate(Id2 = as.numeric(Id2))## Parsed with column specification:
## cols(
## .default = col_double(),
## Id = col_character(),
## Geography = col_character()
## )## See spec(...) for full column specifications.#Read in Freedom Data
freedom <- read_excel("Freedom_In_The_50_States_2018.xlsx", sheet = "Overall")
#Combine
census_edu_2015 <- left_join(census_2015, edu_data_2015, by = c("CensusTract" = "Id2"))
full_2015 <- full_join(census_edu_2015, freedom %>% filter(Year == 2015), by = "State")
#For the education data, we will need to divide the count by the total
full_2015_eduProp <- full_2015 %>% mutate(PropHighSchool = EstHighSchool/EstTotal,
PropBachelors = EstBachelors/EstTotal,
PropDoctorate = EstDoctorate/EstTotal)
full_2015_varOfInterest <- full_2015 %>% select(IncomePerCap,Hispanic:Asian,Professional:Production, SelfEmployed, Unemployment, FiscalPolicy,EconomicFreedom, RegulatoryPolicy, PersonalFreedom, OverallFreedom)“First Attempts” Geo Plots
full_2015_2 <- subset(full_2015, INTPTLONG > -150 & INTPTLONG < -50 & INTPTLAT > 24 & INTPTLAT < 50)
loadPkg("ggplot2")
ggplot(full_2015_2, aes( x= INTPTLONG, y = INTPTLAT, color = Hispanic)) +
geom_point(size=0.005) + scale_color_gradient(low="orange", high="red") +
labs (title="Scatter Plot: Spanish Ethnicity in 2015",
x=" Latitude",
y = "Longitude")
ggplot(full_2015_eduProp, aes( x= Professional, y = IncomePerCap, color = State == "New York", size = TotalPop)) +
geom_point() +
labs (title="Scatter Plot: IncomePerCap and Work type: Professional",
x="Professional (% of pop in census tract)",
y = "IncomePerCap")
# Population Histogram and QQ
loadPkg("ggplot2")
# NA's
dat0<-full_2015[!is.na(full_2015$IncomePerCap),]
sum(is.na(dat0$IncomePerCap))## [1] 0#removed outliers
dat1 <- outlierKD2(dat0, IncomePerCap, TRUE)
## Outliers identified: 3589
## Propotion (%) of outliers: 5.2
## Mean of the outliers: 74140.26
## Mean without removing outliers: 28491.23
## Mean if we remove outliers: 26139.73
## Outliers successfully removed#Histogram of dependent variable
hist(dat1$TotalPop,
col = "#a8c2fb",
main = "population Count",
xlab = "US Dollars")
#ggplot
qqnorm(dat1$TotalPop, main="Q-Q plot of TotalPop")
qqline(dat1$TotalPop)
#Mean after removing NA's
format(mean(dat1$TotalPop, na.rm = TRUE))## [1] "4369.254"#Summary after removing NA's
summary(dat1$TotalPop)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5 2929 4086 4369 5460 53812#standard deviation
sd(dat1$TotalPop, na.rm = TRUE)## [1] 2095.011
Income Histogram and QQ
#Histogram of dependent variable
hist(dat0$IncomePerCap,
col = "#a8c2fb",
main = "Income per Capita",
xlab = "US Dollars")
#ggplot
qqnorm(dat0$IncomePerCap, main="Q-Q plot of IncomePerCap")
qqline(dat0$IncomePerCap)
#Mean after removing NA's
format(mean(dat0$IncomePerCap, na.rm = TRUE))## [1] "28491.23"#Summary after removing NA's
summary(dat0$IncomePerCap)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 128 19123 25344 28491 33894 254204#standard deviation
sd(dat0$IncomePerCap, na.rm = TRUE)## [1] 15047.07#Histogram of dependent variable
hist(dat1$IncomePerCap,
col = "#a8c2fb",
main = "Income per Capita",
xlab = "US Dollars")
#ggplot
qqnorm(dat1$IncomePerCap, main="Q-Q plot of IncomePerCap")
qqline(dat1$IncomePerCap)
#Mean after removing NA's
format(mean(dat1$IncomePerCap, na.rm = TRUE))## [1] "26139.73"#Summary after removing NA's
summary(dat1$IncomePerCap)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 128 18776 24730 26140 32247 56040 3589#standard deviation
sd(dat1$IncomePerCap, na.rm = TRUE)## [1] 10274.98EDA By Group
loadPkg("dvmisc")
#Group by Economic Freedom
GroupEF <- quant_groups(dat1$EconomicFreedom, 4)## Observations per group: 18417, 17756, 22780, 13244. 1064 missing.str(GroupEF)## Factor w/ 4 levels "[-0.827,-0.256]",..: 3 3 3 3 3 3 3 3 3 3 ...#Boxplot of IncomePerCap by Economic Freedom
plot(IncomePerCap ~ GroupEF, data=dat1, main="IncomePerCap and GroupEF", col=c("#ffd18b","#ffb76d","#ffa264", "#ff875f") )
summary(dat1$EconomicFreedom)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.8272 -0.2556 0.0152 -0.0866 0.1376 0.3550 1064#Histogram of Economic Freedom
hist(dat1$EconomicFreedom,
col = "#ffd18b",
main = "Economic Freedom")
#qqplot of Econoimc Freedom
qqnorm(dat1$EconomicFreedom, main="Q-Q plot of Economic Freedom")
qqline(dat1$EconomicFreedom)
#Group by PersonalFreedom
GroupPF <- quant_groups(dat1$PersonalFreedom, 4)## Observations per group: 18061, 18704, 20888, 14544. 1064 missing.str(GroupPF)## Factor w/ 4 levels "[-0.0444,0.0135]",..: 1 1 1 1 1 1 1 1 1 1 ...#Boxplot of IncomePerCap by Personal Freedom
plot(IncomePerCap ~ GroupPF, data=dat1, main="IncomePerCap and GroupPF", col=c("#BFBFFF","#A3A3FF", "#7879FF", "#4949FF") )
summary(dat1$PersonalFreedom)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.0444 0.0135 0.0803 0.0641 0.1064 0.2450 1064#Histogram of Personal Freedom
hist(dat1$PersonalFreedom,
col = "#A3A3FF",
main = "Personal Freedom")
#qqplot of Personal Freedom
qqnorm(dat1$PersonalFreedom, main="Q-Q plot of Personal Freedom")
qqline(dat1$PersonalFreedom)
#Group by Regulatory Policy
GroupRP <- quant_groups(dat1$RegulatoryPolicy, 4)## Observations per group: 19180, 17739, 17836, 17442. 1064 missing.str(GroupRP)## Factor w/ 4 levels "[-0.457,-0.223]",..: 3 3 3 3 3 3 3 3 3 3 ...#Boxplot of IncomePerCap by Regulatory Policy
plot(IncomePerCap ~ GroupRP, data=dat1, main="IncomePerCap and GroupRP", col=c("#FFDF01","#FED901", "#FFCF00", "#FEC300") )
summary(dat1$RegulatoryPolicy)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.4569 -0.2228 -0.0737 -0.1511 -0.0322 0.0715 1064#Histogram of Regulatory Policy
hist(dat1$RegulatoryPolicy,
col = "#FFDF01",
main = "Regulatory Policy ")
#qqplot of Regulatory Policy
qqnorm(dat1$RegulatoryPolicy, main="Q-Q plot of Regulatory Policy")
qqline(dat1$RegulatoryPolicy)
#Group by Fiscal Policy
GroupFP <- quant_groups(dat1$FiscalPolicy, 4)## Observations per group: 18164, 19958, 16580, 17495. 1064 missing.str(GroupFP)## Factor w/ 4 levels "[-0.37,-0.0602]",..: 3 3 3 3 3 3 3 3 3 3 ...#Boxplot of IncomePerCap by Fiscal Policy
plot(IncomePerCap ~ GroupFP, data=dat1, main="IncomePerCap and GroupFP", col=c("#7EFFD4","#70EBBA", "#64D8A7", "#58BD95") )
summary(dat1$FiscalPolicy)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.3702 -0.0602 0.0634 0.0646 0.1767 0.4024 1064#Histogram of Fiscal Policy
hist(dat1$FiscalPolicy,
col = "#7EFFD4",
main = "Fiscal Policy")
#qqplot of Fiscal Policy
qqnorm(dat1$FiscalPolicy, main="Q-Q plot of Fiscal Policy")
qqline(dat1$FiscalPolicy)
#Group by Overall Freedom
GroupOF <- quant_groups(dat1$OverallFreedom, 4)## Observations per group: 18440, 17906, 18395, 17456. 1064 missing.str(GroupOF)## Factor w/ 4 levels "[-0.814,-0.102]",..: 2 2 2 2 2 2 2 2 2 2 ...#Boxplot of IncomePerCap by Overall Freedom
plot(IncomePerCap ~ GroupOF, data=dat1, main="IncomePerCap and GroupOF", col=c("#C0C6CB","#98A5C0", "#7688BB", "#536CB5") )
summary(dat1$OverallFreedom)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.8136 -0.1021 0.0652 -0.0224 0.1637 0.4614 1064#Histogram of Overall Freedom
hist(dat1$OverallFreedom,
col = "#A3A3FF",
main = "Overall Freedom")
#qqplot of Overall Freedom
qqnorm(dat1$OverallFreedom, main="Q-Q plot of Overall Freedom")
qqline(dat1$OverallFreedom)
#Group by Unemployment
GroupUnemployment <- quant_groups(dat1$Unemployment, 4)## Observations per group: 18765, 18330, 18095, 17970. 101 missing.str(GroupUnemployment)## Factor w/ 4 levels "[0,5.1]","(5.1,7.7]",..: 2 4 2 3 1 3 3 3 3 2 ...#Boxplot of IncomePerCap by Unemployment
plot(IncomePerCap ~ GroupUnemployment, data=dat1, main="IncomePerCap and GroupUnemployment", col=c("#FFF57B","#FFE469", "#FECC51", "#FCB033") )
summary(dat1$Unemployment)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 5.100 7.700 9.028 11.400 100.000 101#Histogram of Unemployment
hist(dat1$Unemployment,
col = "#A3A3FF",
)
#qqplot of Unemployment Unemployment
qqnorm(dat1$Unemployment, main="Q-Q plot of Unemployment")
qqline(dat1$Unemployment)
#Group by Professional
GroupProfessional <- quant_groups(dat1$Professional, 4)## Observations per group: 18379, 18323, 18173, 18281. 105 missing.str(GroupProfessional)## Factor w/ 4 levels "[0,24.1]","(24.1,32.6]",..: 3 1 2 2 4 2 1 3 2 2 ...#Boxplot of IncomePerCap by Professional
plot(IncomePerCap ~ GroupProfessional, data=dat1, main="IncomePerCap and GroupProfessional", col=c("#F3D8F2","#E6B2E4", "#D88DD5", "#CA68C7") )
summary(dat1$Professional)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 24.1 32.6 34.8 43.8 100.0 105#Histogram of Professional
hist(dat1$Professional,
col = "#D88DD5",
main = "Professional")
#qqplot of Professional
qqnorm(dat1$Professional, main="Q-Q plot of Professional")
qqline(dat1$Professional)
#Group by Office
GroupOffice <- quant_groups(dat1$Office, 4)## Observations per group: 18303, 18828, 17766, 18259. 105 missing.str(GroupOffice)## Factor w/ 4 levels "[0,20.1]","(20.1,23.8]",..: 2 2 2 3 1 4 3 4 3 1 ...#Boxplot of IncomePerCap by Office
plot(IncomePerCap ~ GroupOffice, data=dat1, main="IncomePerCap and GroupOffice", col=c("#D1FFD5","#B4FFB2", "#98FF98", "#79F58A") )
summary(dat1$Office)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 20.10 23.80 23.95 27.50 100.00 105#Histogram of Office
hist(dat1$Office,
col = "#FECC51",
main = "Office")
#qqplot of Office
qqnorm(dat1$Office, main="Q-Q plot of Office")
qqline(dat1$Office)
#Group by Service
GroupService <- quant_groups(dat1$Service, 4)## Observations per group: 18672, 18068, 18275, 18141. 105 missing.str(GroupService)## Factor w/ 4 levels "[0,13.5]","(13.5,17.9]",..: 2 4 4 3 2 2 4 1 2 2 ...#Boxplot of IncomePerCap by Service
plot(IncomePerCap ~ GroupService, data=dat1, main="IncomePerCap and GroupService", col=c("#E0BCBF","#D8ABB1", "#CF989F", "#C0838C") )
summary(dat1$Service)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 13.5 17.9 19.1 23.6 100.0 105#Histogram of Service
hist(dat1$Service,
col = "#79F58A",
main = "Service")
#qqplot of Service
qqnorm(dat1$Service, main="Q-Q plot of Service")
qqline(dat1$Service)
#Group by Construction
GroupConstruction <- quant_groups(dat1$Construction, 4)## Observations per group: 18588, 18209, 18113, 18246. 105 missing.str(GroupConstruction)## Factor w/ 4 levels "[0,5]","(5,8.4]",..: 3 3 3 3 1 2 3 2 2 3 ...#Boxplot of IncomePerCap by Construction
plot(IncomePerCap ~ GroupConstruction, data=dat1, main="IncomePerCap and GroupConstruction", col=c("#A9AB98","#949180", "#7D7968", "#5E594F") )
summary(dat1$Construction)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 5.000 8.400 9.295 12.500 100.000 105#Histogram of Construciton
hist(dat1$Construction,
col = "#C0838C",
main = "Construction")
#qqplot of Construction
qqnorm(dat1$Construction, main="Q-Q plot of Construction")
qqline(dat1$Construction)
#Group by Production
GroupProduction <- quant_groups(dat1$Production, 4)## Observations per group: 18566, 18226, 18085, 18279. 105 missing.str(GroupProduction)## Factor w/ 4 levels "[0,7.1]","(7.1,11.8]",..: 3 4 3 3 3 3 3 1 3 4 ...#Boxplot of IncomePerCap by Production
plot(IncomePerCap ~ GroupProduction, data=dat1, main="IncomePerCap and Production", col=c("#F3D8F2","#E6B2E4", "#D88DD5", "#CA68C7") )
summary(dat1$Production)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 7.10 11.80 12.86 17.40 100.00 105#Histogram of Producrtion
hist(dat1$Production,
col = "#E6B2E4",
main = "Production")
#qqplot of Production
qqnorm(dat1$Production, main="Q-Q plot of Production")
qqline(dat1$Production)
#Group by Self-Employed
GroupProduction <- quant_groups(dat1$SelfEmployed, 4)## Observations per group: 19155, 17839, 18195, 17967. 105 missing.str(GroupProduction)## Factor w/ 4 levels "[0,3.6]","(3.6,5.5]",..: 2 3 4 1 2 3 1 3 2 3 ...#Boxplot of IncomePerCap by Selfemployed
plot(IncomePerCap ~ GroupProduction, data=dat1, main="IncomePerCap and Selfemployed", col=c("#F3D8F2","#E6B2E4", "#D88DD5", "#CA68C7") )
summary(dat1$SelfEmployed)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 3.600 5.500 6.227 8.100 100.000 105#Histogram of Self Employed
hist(dat1$SelfEmployed,
col = "#E6B2E4",
main = "Production")
#qqplot of Self employed
qqnorm(dat1$SelfEmployed, main="Q-Q plot of Selfemployed")
qqline(dat1$SelfEmployed)
#Group by Black
GroupBlack <- quant_groups(dat1$Black, 4)## Observations per group: 18430, 18252, 18303, 18276. 0 missing.str(GroupBlack)## Factor w/ 4 levels "[0,0.7]","(0.7,3.7]",..: 3 4 4 2 4 3 4 3 3 3 ...#Boxplot of IncomePerCap by Black Population
plot(IncomePerCap ~ GroupBlack, data=dat1, main="IncomePerCap and GroupBlack", col=c("#d0dfff","#a8c2fb", "#86abf9", "#6893ee") )
summary(dat1$Black)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.70 3.70 13.27 14.40 100.00#Histogram of Black
hist(dat1$Black,
col = "#d0dfff",
main = "Black")
#qqplot of Black
qqnorm(dat1$Black, main="Q-Q plot of Black")
qqline(dat1$Black)
#Group by Hispanic
GroupHispanic <- quant_groups(dat1$Hispanic, 4)## Observations per group: 18472, 18246, 18233, 18310. 0 missing.str(GroupHispanic)## Factor w/ 4 levels "[0,2.4]","(2.4,7]",..: 1 1 1 3 1 3 2 1 1 1 ...#Boxplot of IncomePerCap by Hispanic Population
plot(IncomePerCap ~ GroupHispanic, data=dat1, main="IncomePerCap and GroupHispanic", col=c("#F6BDC0","#F1959B", "#F07470", "#EA4C46") )
summary(dat1$Hispanic)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.40 7.00 16.86 20.40 100.00#Histogram of Hispanic
hist(dat1$Hispanic,
col = "#EA4C46",
main = "Hispanic")
#qqplot of Hispanic
qqnorm(dat1$Hispanic, main="Q-Q plot of Hspanic")
qqline(dat1$Hispanic)
#Group by Asians
GroupAsian <- quant_groups(dat1$Asian, 3)## Observations per group: 26078, 23036, 24147. 0 missing.str(GroupAsian)## Factor w/ 3 levels "[0,0.5]","(0.5,3.2]",..: 2 2 2 1 3 1 1 1 1 1 ...#Boxplot of IncomePerCap by Asian Population
plot(IncomePerCap ~ GroupAsian, data=dat1, main="IncomePerCap and GroupAsian", col=c("#FFE6C8","#FFCCBE", "#EFABA0", "#D6806F") )
summary(dat1$Asian)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.20 1.40 4.59 4.80 91.30#Histogram of Asian
hist(dat1$Asian,
col = "#FFE6C8",
main = "Asian")
#qqplot of Asian
qqnorm(dat1$Asian, main="Q-Q plot of Asian")
qqline(dat1$Asian)
#Group by White
GroupWhite <- quant_groups(dat1$White, 4)## Observations per group: 18351, 18331, 18285, 18294. 0 missing.str(GroupWhite)## Factor w/ 4 levels "[0,39.4]","(39.4,71.4]",..: 3 2 3 3 2 3 3 3 4 3 ...#Boxplot of IncomePerCap by White Population
plot(IncomePerCap ~ GroupWhite, data=dat1, main="IncomePerCap and GroupWhite", col=c("#F5F8FA", "#EFF2F4", "#EAEDEF", "#E0E3E5") )
summary(dat1$White)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 39.40 71.40 62.03 88.30 100.00#Histogram of White
hist(dat1$White,
col = "#F5F8FA",
main = "White")
#qqplot of White
qqnorm(dat1$White, main="Q-Q plot of White")
qqline(dat1$White)
Correlations
full_2015_cor <- cor(full_2015_varOfInterest, use = "complete.obs")
loadPkg("corrplot")## corrplot 0.84 loadedcorrplot(full_2015_cor)
Freedom_2015 <- subset(full_2015_varOfInterest, select = c(IncomePerCap, OverallFreedom, EconomicFreedom, PersonalFreedom, RegulatoryPolicy, FiscalPolicy))
Work_2015 <- subset(full_2015_varOfInterest, select = c(IncomePerCap, Professional, Service, Office, Construction, Production, SelfEmployed, Unemployment))
Ethnic_2015 <- subset(full_2015_varOfInterest, select = c(IncomePerCap, Hispanic, White, Black, Native, Asian))Freedom_2015_cor <- cor(Freedom_2015, use = "complete.obs")
#corrplot(Freedom_2015_cor)
corrplot.mixed(Freedom_2015_cor, tl.pos = "lt")
Work_2015_cor <- cor(Work_2015, use = "complete.obs")
#corrplot(Work_2015_cor)
corrplot.mixed(Work_2015_cor, tl.pos = "lt")
Ethnic_2015_cor <- cor(Ethnic_2015, use = "complete.obs")
#corrplot(Ethnic_2015_cor)
corrplot.mixed(Ethnic_2015_cor, tl.pos = "lt")
Scatter PLots
ggplot(dat1, aes( x= Professional, y = IncomePerCap, color = White)) +
geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~Professional) +
labs (title="Scatter Plot of IncomePerCap and Professional",
x="Professional",
y = "IncomePerCap")
ggplot(dat1, aes( x= Unemployment, y = IncomePerCap, color = White)) +
geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~Unemployment) +
labs (title="Scatter Plot of IncomePerCap and Unemployment",
x="Unemployment",
y = "IncomePerCap")
ggplot(dat1, aes( x= Service, y = IncomePerCap, color = White)) +
geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~Service) +
labs (title="Scatter Plot of IncomePerCap and Service",
x="Service",
y = "IncomePerCap")
ggplot(dat1, aes( x= EconomicFreedom, y = IncomePerCap, color = White)) +
geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~EconomicFreedom) +
labs (title="Scatter Plot of IncomePerCap and EconomicFreedom",
x="EconomicFreedom",
y = "IncomePerCap")
ANOVA
#' @title Returns the max value of each row of a data.frame or matrix
#'
#' @description
#' Returns maximum value of each row of a data.frame or matrix.
#' @param df Data.frame or matrix, required.
#' @param na.rm Logical value, optional, TRUE by default. Defines whether NA values should be removed first. Otherwise result will be NA when any NA is in the given vector.
#' @return Returns a vector of numbers of length equal to number of rows in df.
#' @template maxmin
#' @export
rowMaxs <- function(df, na.rm=TRUE) {
if (is.matrix(df)) {df <- data.frame(df, stringsAsFactors=FALSE, drop = FALSE)}
valid.cols <- sapply(df, function(x) { is.numeric(x) || is.logical(x) || is.character(x)})
stopifnot(any(valid.cols))
# or could just return NA?:
# if (!any(valid.cols)) {return(NA)}
if (any(!valid.cols) ) {warning('using only numeric (double or integer) or logical or character columns -- ignoring other columns ')}
result <- do.call(pmax, c(df[ , valid.cols, drop = FALSE], na.rm=na.rm))
result[nononmissing <- rowSums(!is.na(df[ , valid.cols, drop = FALSE]))==0] <- -Inf
if (any(nononmissing)) {warning('where no non-missing arguments, returning -Inf')}
return(result)
# df = data.frame of numeric values, i.e. a list of vectors passed to pmax
# Value returned is vector, each element is max of a row of df
}
full_2015_varOfInterest1 <- full_2015_varOfInterest %>% mutate(EthnicMajority = case_when(
White >= 50 ~ "White",
Hispanic >= 50 ~ "Hispanic",
Black >= 50 ~ "Black",
Native >= 50 ~ "Native",
Asian >= 50 ~ "Asian",
TRUE ~ "Diverse"))
EthnicityMaxes <- rowMaxs(full_2015_varOfInterest1 %>% select(White, Hispanic, Black, Native, Asian))
WorkMaxes <- rowMaxs(full_2015_varOfInterest1 %>% select(Professional:Unemployment))
full_2015_varOfInterest2 <- cbind(full_2015_varOfInterest1, EthnicityMaxes, WorkMaxes)
full_2015_varOfInterest_withMajorityEthnicity <- full_2015_varOfInterest2 %>% mutate(EthnicPlurality = case_when(
White == EthnicityMaxes ~ "White",
Hispanic == EthnicityMaxes ~ "Hispanic",
Black == EthnicityMaxes ~ "Black",
Native == EthnicityMaxes ~ "Native",
Asian == EthnicityMaxes ~ "Asian",
TRUE ~ "Error"),
WorkPlurality = case_when(
Professional == WorkMaxes ~ "Professional",
Service == WorkMaxes ~ "Service",
Office == WorkMaxes ~ "Office",
Construction == WorkMaxes ~ "Construction",
Production == WorkMaxes ~ "Production",
SelfEmployed == WorkMaxes ~ "SelfEmployed",
Unemployment == WorkMaxes ~ "Unemployment",
TRUE ~ "Error"))
full_2015_varOfInterest_withMajorityEthnicity$EthnicMajority <-
as.factor(full_2015_varOfInterest_withMajorityEthnicity$EthnicMajority)
full_2015_varOfInterest_withMajorityEthnicity$EthnicPlurality <-
as.factor(full_2015_varOfInterest_withMajorityEthnicity$EthnicPlurality)
full_2015_varOfInterest_withMajorityEthnicity$WorkPlurality <-
as.factor(full_2015_varOfInterest_withMajorityEthnicity$WorkPlurality)anova_dat <- full_2015_varOfInterest_withMajorityEthnicity %>% filter(EthnicPlurality!= "Error",
WorkPlurality!= "Error")Ethnic_2015_anova = aov(IncomePerCap ~ EthnicPlurality, data = anova_dat)
Ethnic_2015_anova## Call:
## aov(formula = IncomePerCap ~ EthnicPlurality, data = anova_dat)
##
## Terms:
## EthnicPlurality Residuals
## Sum of Squares 2.590039e+12 1.394545e+13
## Deg. of Freedom 4 73155
##
## Residual standard error: 13806.84
## Estimated effects may be unbalanced
## 39 observations deleted due to missingnessnames(Ethnic_2015_anova)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "na.action" "contrasts" "xlevels" "call"
## [13] "terms" "model"summary(Ethnic_2015_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## EthnicPlurality 4 2.590e+12 6.475e+11 3397 <2e-16 ***
## Residuals 73155 1.395e+13 1.906e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 39 observations deleted due to missingnessggplot(anova_dat, aes(x=EthnicPlurality, y=IncomePerCap)) +
geom_boxplot(outlier.shape=8, outlier.size=4) +
labs(title="Income/Capita with Different Ethnic Majority", x="Ethnic Majority", y = "Income per Cap")
TukeyHSD(Ethnic_2015_anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = IncomePerCap ~ EthnicPlurality, data = anova_dat)
##
## $EthnicPlurality
## diff lwr upr p adj
## Black-Asian -12780.4809 -13907.7147 -11653.24708 0.0000000
## Hispanic-Asian -13705.6676 -14812.6303 -12598.70484 0.0000000
## Native-Asian -16464.3603 -19310.6485 -13618.07205 0.0000000
## White-Asian 647.5161 -403.9082 1698.94040 0.4464321
## Hispanic-Black -925.1867 -1505.7776 -344.59577 0.0001346
## Native-Black -3683.8794 -6369.5964 -998.16235 0.0017134
## White-Black 13427.9970 12961.9366 13894.05745 0.0000000
## Native-Hispanic -2758.6927 -5435.9649 -81.42052 0.0396499
## White-Hispanic 14353.1837 13938.5480 14767.81938 0.0000000
## White-Native 17111.8764 14457.0858 19766.66697 0.0000000Work_2015_anova = aov(IncomePerCap ~ WorkPlurality, data = anova_dat)
Work_2015_anova## Call:
## aov(formula = IncomePerCap ~ WorkPlurality, data = anova_dat)
##
## Terms:
## WorkPlurality Residuals
## Sum of Squares 4.242548e+12 1.229294e+13
## Deg. of Freedom 6 73153
##
## Residual standard error: 12963.19
## Estimated effects may be unbalanced
## 39 observations deleted due to missingnessnames(Work_2015_anova)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "na.action" "contrasts" "xlevels" "call"
## [13] "terms" "model"summary(Work_2015_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## WorkPlurality 6 4.243e+12 7.071e+11 4208 <2e-16 ***
## Residuals 73153 1.229e+13 1.680e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 39 observations deleted due to missingnessggplot(anova_dat, aes(x=WorkPlurality, y=IncomePerCap)) +
geom_boxplot(outlier.shape=8, outlier.size=4) +
labs(title="Income/Capita with Different Work Majority", x="Work Majority", y = "Income per Cap")
TukeyHSD(Work_2015_anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = IncomePerCap ~ WorkPlurality, data = anova_dat)
##
## $WorkPlurality
## diff lwr upr p adj
## Office-Construction 4152.88357 2905.526990 5400.2402 0.0000000
## Production-Construction 1332.84270 -5.891578 2671.5770 0.0518977
## Professional-Construction 17590.26271 16385.932493 18794.5929 0.0000000
## SelfEmployed-Construction 11446.99976 349.799288 22544.2002 0.0380499
## Service-Construction 73.21197 -1179.886893 1326.3108 0.9999978
## Unemployment-Construction -5476.10244 -7681.711665 -3270.4932 0.0000000
## Production-Office -2820.04087 -3533.458145 -2106.6236 0.0000000
## Professional-Office 13437.37913 13028.519743 13846.2385 0.0000000
## SelfEmployed-Office 7294.11618 -3745.114449 18333.3468 0.4482729
## Service-Office -4079.67161 -4615.406142 -3543.9371 0.0000000
## Unemployment-Office -9628.98602 -11521.462377 -7736.5097 0.0000000
## Professional-Production 16257.42001 15622.221597 16892.6184 0.0000000
## SelfEmployed-Production 10114.15706 -935.771630 21164.0857 0.0985074
## Service-Production -1259.63073 -1983.041067 -536.2204 0.0000059
## Unemployment-Production -6808.94514 -8762.858599 -4855.0317 0.0000000
## SelfEmployed-Professional -6143.26295 -17177.714716 4891.1888 0.6552473
## Service-Professional -17517.05074 -17943.107437 -17090.9940 0.0000000
## Unemployment-Professional -23066.36515 -24930.763066 -21201.9672 0.0000000
## Service-SelfEmployed -11373.78779 -22413.668735 -333.9068 0.0384811
## Unemployment-SelfEmployed -16923.10220 -28111.195272 -5735.0091 0.0001665
## Unemployment-Service -5549.31441 -7445.580501 -3653.0483 0.0000000Chi Squared Tests
loadPkg('sjPlot')## Learn more about sjPlot with 'browseVignettes("sjPlot")'.loadPkg('visualize')
chi_table_dat <- full_2015_varOfInterest_withMajorityEthnicity %>%
select(c(WorkPlurality, EthnicPlurality)) %>% filter(WorkPlurality != "Error")
chi_table_dat %>%
sjtab(fun = "xtab", var.labels=c("Work", "Ethnicity"),
show.summary=T, show.exp=T, show.legend=T)| Work | Ethnicity | Total | |||||
|---|---|---|---|---|---|---|---|
| Asian | Black | Error | Hispanic | Native | White | ||
| Construction |
0 19 |
31 104 |
0 0 |
639 137 |
1 3 |
359 767 |
1030 1030 |
| Error |
0 0 |
0 0 |
0 0 |
0 0 |
0 0 |
0 0 |
0 0 |
| Office |
128 193 |
1564 1087 |
0 0 |
2199 1426 |
17 30 |
6816 7989 |
10724 10724 |
| Production |
19 70 |
480 398 |
0 0 |
993 522 |
6 11 |
2426 2923 |
3924 3924 |
| Professional |
927 851 |
2346 4804 |
0 0 |
2510 6299 |
114 131 |
41484 35296 |
47381 47381 |
| SelfEmployed |
0 0 |
0 1 |
0 0 |
2 2 |
0 0 |
11 10 |
13 13 |
| Service |
241 174 |
2744 984 |
0 0 |
3274 1290 |
49 27 |
3394 7227 |
9702 9702 |
| Unemployment |
0 8 |
257 43 |
0 0 |
114 56 |
15 1 |
39 317 |
425 425 |
| Total |
1315 1315 |
7422 7422 |
0 0 |
9731 9731 |
202 202 |
54529 54529 |
73199 73199 |
χ2=NaN · df=35 · Cramer’s V=NaN · Fisher’s p=0.000 |
observed values
expected values
chi.test <- chisq.test(full_2015_varOfInterest_withMajorityEthnicity$EthnicPlurality, full_2015_varOfInterest_withMajorityEthnicity$WorkPlurality)#Exhaustive Search
loadPkg('leaps')
#This is essentially best fit
reg.best10 <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17) ## Reordering variables and trying again:plot(reg.best10, scale = "adjr2", main = "Adjusted R^2")
plot(reg.best10, scale = "r2", main = "R^2")
# In the "leaps" package, we can use scale=c("bic","Cp","adjr2","r2")
plot(reg.best10, scale = "bic", main = "BIC")
plot(reg.best10, scale = "Cp", main = "Cp")
reg.summary<-summary(reg.best10)
names(reg.best10)## [1] "np" "nrbar" "d" "rbar" "thetab" "first"
## [7] "last" "vorder" "tol" "rss" "bound" "nvmax"
## [13] "ress" "ir" "nbest" "lopt" "il" "ier"
## [19] "xnames" "method" "force.in" "force.out" "sserr" "intercept"
## [25] "lindep" "reorder" "nullrss" "nn" "call"plot(reg.summary$rsq,xlab = "# of Varibales", ylab = "Rsquare", type = "l")
plot(reg.summary$rss,xlab = "# of Varibales", ylab = "RSS", type = "l")
#Adjusted R2
which.max(reg.summary$adjr2)## [1] 14plot(reg.summary$adjr2,xlab = '# of Varibales', ylab = 'Adjusted Rsq', type = "l")
points(12,reg.summary$adjr2[12],col = "red",cex = 2, pch =20)
#CP
which.min(reg.summary$cp)## [1] 13plot(reg.summary$cp,xlab = '# of Varibales', ylab = 'CP', type = "l")
points(12,reg.summary$cp[12],col = "red",cex = 2, pch =20)
#BIC
which.min(reg.summary$bic)## [1] 12plot(reg.summary$bic,xlab = '# of Varibales', ylab = 'BIC', type = "l")
points(10,reg.summary$bic[12],col = "red",cex = 2, pch =20)
Forward Selection
reg.forward10 <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17, nbest = 2, method = "forward")## Reordering variables and trying again:plot(reg.forward10, scale = "adjr2", main = "Adjusted R^2")
plot(reg.forward10, scale = "bic", main = "BIC")
plot(reg.forward10, scale = "Cp", main = "Cp")
#summary(reg.forward10)Backward Selection
Now backwards (nvmax=10 and nbest=2)
reg.back10 <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17, nbest = 2, method = "backward")## Reordering variables and trying again:plot(reg.back10, scale = "adjr2", main = "Adjusted R^2")
plot(reg.back10, scale = "bic", main = "BIC")
plot(reg.back10, scale = "Cp", main = "Cp")
summary(reg.back10)## Subset selection object
## Call: regsubsets.formula(IncomePerCap ~ ., data = full_2015_varOfInterest,
## nvmax = 17, nbest = 2, method = "backward")
## 17 Variables (and intercept)
## Forced in Forced out
## Hispanic FALSE FALSE
## White FALSE FALSE
## Black FALSE FALSE
## Native FALSE FALSE
## Asian FALSE FALSE
## Professional FALSE FALSE
## Service FALSE FALSE
## Office FALSE FALSE
## Construction FALSE FALSE
## Production FALSE FALSE
## SelfEmployed FALSE FALSE
## Unemployment FALSE FALSE
## FiscalPolicy FALSE FALSE
## EconomicFreedom FALSE FALSE
## PersonalFreedom FALSE FALSE
## RegulatoryPolicy FALSE FALSE
## OverallFreedom FALSE FALSE
## 2 subsets of each size up to 15
## Selection Algorithm: backward
## Hispanic White Black Native Asian Professional Service Office
## 1 ( 1 ) " " " " " " " " " " "*" " " " "
## 1 ( 2 ) " " "*" " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " "*" " " " "
## 2 ( 2 ) " " "*" "*" " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " "*" " " " "
## 3 ( 2 ) " " "*" " " " " " " "*" " " " "
## 4 ( 1 ) " " "*" " " " " " " "*" " " " "
## 4 ( 2 ) " " "*" " " " " " " "*" " " " "
## 5 ( 1 ) " " "*" " " " " " " "*" " " " "
## 5 ( 2 ) " " "*" " " " " " " "*" "*" " "
## 6 ( 1 ) " " "*" " " " " " " "*" " " " "
## 6 ( 2 ) " " "*" " " " " " " "*" "*" " "
## 7 ( 1 ) " " "*" " " " " " " "*" "*" " "
## 7 ( 2 ) " " "*" " " " " " " "*" "*" " "
## 8 ( 1 ) " " "*" " " " " " " "*" "*" " "
## 8 ( 2 ) " " "*" " " " " " " "*" "*" " "
## 9 ( 1 ) " " "*" " " " " " " "*" "*" " "
## 9 ( 2 ) " " "*" "*" " " " " "*" "*" " "
## 10 ( 1 ) " " "*" "*" " " " " "*" "*" " "
## 10 ( 2 ) "*" "*" "*" " " " " "*" "*" " "
## 11 ( 1 ) "*" "*" "*" " " " " "*" "*" " "
## 11 ( 2 ) "*" "*" "*" "*" " " "*" "*" " "
## 12 ( 1 ) "*" "*" "*" "*" " " "*" "*" " "
## 12 ( 2 ) "*" "*" "*" "*" "*" "*" "*" " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 13 ( 2 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 2 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## Construction Production SelfEmployed Unemployment FiscalPolicy
## 1 ( 1 ) " " " " " " " " " "
## 1 ( 2 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " "
## 2 ( 2 ) " " " " " " " " " "
## 3 ( 1 ) " " " " "*" " " " "
## 3 ( 2 ) " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " "*"
## 4 ( 2 ) " " " " "*" "*" " "
## 5 ( 1 ) " " " " " " "*" "*"
## 5 ( 2 ) " " " " "*" "*" " "
## 6 ( 1 ) " " " " "*" "*" "*"
## 6 ( 2 ) " " " " "*" "*" "*"
## 7 ( 1 ) " " " " "*" "*" "*"
## 7 ( 2 ) "*" "*" "*" "*" " "
## 8 ( 1 ) "*" " " "*" "*" "*"
## 8 ( 2 ) "*" "*" "*" "*" "*"
## 9 ( 1 ) "*" "*" "*" "*" "*"
## 9 ( 2 ) "*" "*" "*" "*" "*"
## 10 ( 1 ) "*" "*" "*" "*" "*"
## 10 ( 2 ) "*" "*" "*" "*" "*"
## 11 ( 1 ) "*" "*" "*" "*" "*"
## 11 ( 2 ) "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" "*"
## 12 ( 2 ) "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*" "*"
## 13 ( 2 ) "*" "*" "*" "*" "*"
## 14 ( 1 ) "*" "*" "*" "*" "*"
## 14 ( 2 ) "*" "*" "*" "*" "*"
## 15 ( 1 ) "*" "*" "*" "*" "*"
## EconomicFreedom RegulatoryPolicy PersonalFreedom OverallFreedom
## 1 ( 1 ) " " " " " " " "
## 1 ( 2 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 2 ( 2 ) " " " " " " " "
## 3 ( 1 ) " " " " " " " "
## 3 ( 2 ) "*" " " " " " "
## 4 ( 1 ) "*" " " " " " "
## 4 ( 2 ) " " " " " " " "
## 5 ( 1 ) "*" " " " " " "
## 5 ( 2 ) " " " " " " " "
## 6 ( 1 ) "*" " " " " " "
## 6 ( 2 ) " " " " " " " "
## 7 ( 1 ) "*" " " " " " "
## 7 ( 2 ) " " " " " " " "
## 8 ( 1 ) "*" " " " " " "
## 8 ( 2 ) " " " " " " " "
## 9 ( 1 ) "*" " " " " " "
## 9 ( 2 ) " " " " " " " "
## 10 ( 1 ) "*" " " " " " "
## 10 ( 2 ) " " " " " " " "
## 11 ( 1 ) "*" " " " " " "
## 11 ( 2 ) " " " " " " " "
## 12 ( 1 ) "*" " " " " " "
## 12 ( 2 ) " " " " " " " "
## 13 ( 1 ) "*" " " " " " "
## 13 ( 2 ) " " " " " " " "
## 14 ( 1 ) "*" " " "*" " "
## 14 ( 2 ) "*" " " " " " "
## 15 ( 1 ) "*" " " "*" " "Sequential Replacement seqrep
reg.seqrep <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17, nbest = 2 , method = "seqrep")## Reordering variables and trying again:plot(reg.seqrep, scale = "adjr2", main = "Adjusted R^2")
plot(reg.seqrep, scale = "bic", main = "BIC")
plot(reg.seqrep, scale = "Cp", main = "Cp")
summary(reg.seqrep)## Subset selection object
## Call: regsubsets.formula(IncomePerCap ~ ., data = full_2015_varOfInterest,
## nvmax = 17, nbest = 2, method = "seqrep")
## 17 Variables (and intercept)
## Forced in Forced out
## Hispanic FALSE FALSE
## White FALSE FALSE
## Black FALSE FALSE
## Native FALSE FALSE
## Asian FALSE FALSE
## Professional FALSE FALSE
## Service FALSE FALSE
## Office FALSE FALSE
## Construction FALSE FALSE
## Production FALSE FALSE
## SelfEmployed FALSE FALSE
## Unemployment FALSE FALSE
## FiscalPolicy FALSE FALSE
## EconomicFreedom FALSE FALSE
## PersonalFreedom FALSE FALSE
## RegulatoryPolicy FALSE FALSE
## OverallFreedom FALSE FALSE
## 2 subsets of each size up to 15
## Selection Algorithm: 'sequential replacement'
## Hispanic White Black Native Asian Professional Service Office
## 1 ( 1 ) " " " " " " " " " " "*" " " " "
## 1 ( 2 ) " " " " " " " " " " " " "*" " "
## 2 ( 1 ) " " " " " " " " " " "*" " " " "
## 2 ( 2 ) " " " " "*" " " " " "*" " " " "
## 3 ( 1 ) " " "*" " " " " " " "*" " " " "
## 3 ( 2 ) " " " " " " " " " " "*" " " " "
## 4 ( 1 ) " " "*" " " " " " " "*" " " " "
## 4 ( 2 ) " " "*" " " " " " " "*" " " " "
## 5 ( 1 ) " " "*" " " " " " " "*" " " " "
## 5 ( 2 ) " " " " " " " " " " "*" "*" " "
## 6 ( 1 ) " " "*" " " " " " " "*" " " "*"
## 6 ( 2 ) " " "*" " " " " " " "*" "*" " "
## 7 ( 1 ) " " "*" " " " " " " "*" " " "*"
## 7 ( 2 ) " " "*" " " " " " " "*" " " "*"
## 8 ( 1 ) " " "*" " " " " " " "*" "*" "*"
## 8 ( 2 ) " " "*" " " " " " " "*" "*" "*"
## 9 ( 1 ) " " "*" " " "*" " " "*" "*" "*"
## 9 ( 2 ) " " "*" " " "*" " " "*" "*" "*"
## 10 ( 1 ) " " "*" " " "*" "*" "*" "*" "*"
## 10 ( 2 ) " " "*" " " "*" "*" "*" "*" "*"
## 11 ( 1 ) " " "*" " " "*" "*" "*" "*" " "
## 11 ( 2 ) " " "*" " " "*" "*" "*" "*" " "
## 12 ( 1 ) "*" "*" "*" "*" " " "*" "*" " "
## 12 ( 2 ) "*" "*" "*" "*" " " "*" "*" " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 13 ( 2 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 2 ) "*" "*" "*" "*" "*" "*" "*" " "
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 15 ( 2 ) "*" "*" "*" "*" "*" "*" "*" "*"
## Construction Production SelfEmployed Unemployment FiscalPolicy
## 1 ( 1 ) " " " " " " " " " "
## 1 ( 2 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " "*" " "
## 2 ( 2 ) " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " "
## 3 ( 2 ) " " " " " " "*" " "
## 4 ( 1 ) " " " " " " "*" " "
## 4 ( 2 ) " " " " "*" " " " "
## 5 ( 1 ) " " " " "*" "*" " "
## 5 ( 2 ) " " " " "*" "*" " "
## 6 ( 1 ) " " " " "*" "*" " "
## 6 ( 2 ) " " " " "*" "*" " "
## 7 ( 1 ) " " " " "*" "*" "*"
## 7 ( 2 ) " " " " "*" "*" " "
## 8 ( 1 ) " " " " "*" "*" " "
## 8 ( 2 ) " " " " "*" "*" "*"
## 9 ( 1 ) " " " " "*" "*" "*"
## 9 ( 2 ) " " " " "*" "*" "*"
## 10 ( 1 ) " " " " "*" "*" " "
## 10 ( 2 ) " " " " "*" "*" "*"
## 11 ( 1 ) "*" "*" "*" "*" " "
## 11 ( 2 ) "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" " "
## 12 ( 2 ) "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*" " "
## 13 ( 2 ) "*" "*" "*" "*" "*"
## 14 ( 1 ) "*" "*" "*" "*" "*"
## 14 ( 2 ) "*" "*" "*" "*" "*"
## 15 ( 1 ) "*" "*" "*" "*" " "
## 15 ( 2 ) "*" "*" "*" "*" " "
## EconomicFreedom RegulatoryPolicy PersonalFreedom OverallFreedom
## 1 ( 1 ) " " " " " " " "
## 1 ( 2 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 2 ( 2 ) " " " " " " " "
## 3 ( 1 ) " " "*" " " " "
## 3 ( 2 ) " " "*" " " " "
## 4 ( 1 ) " " "*" " " " "
## 4 ( 2 ) " " "*" " " " "
## 5 ( 1 ) " " "*" " " " "
## 5 ( 2 ) " " "*" " " " "
## 6 ( 1 ) " " "*" " " " "
## 6 ( 2 ) " " "*" " " " "
## 7 ( 1 ) " " "*" " " " "
## 7 ( 2 ) "*" "*" " " " "
## 8 ( 1 ) "*" "*" " " " "
## 8 ( 2 ) " " "*" " " " "
## 9 ( 1 ) "*" " " " " " "
## 9 ( 2 ) " " "*" " " " "
## 10 ( 1 ) "*" "*" " " " "
## 10 ( 2 ) "*" " " " " " "
## 11 ( 1 ) "*" "*" " " " "
## 11 ( 2 ) " " "*" " " " "
## 12 ( 1 ) "*" "*" " " " "
## 12 ( 2 ) " " "*" " " " "
## 13 ( 1 ) "*" "*" " " " "
## 13 ( 2 ) " " "*" " " " "
## 14 ( 1 ) " " "*" " " "*"
## 14 ( 2 ) " " " " "*" "*"
## 15 ( 1 ) "*" "*" "*" " "
## 15 ( 2 ) "*" "*" " " "*"